Flow-induced channelization in a porous medium 
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We propose a theory for erosional channelization induced by fluid flow in a saturated granular 
porous medium. When the local fluid flow-induced stress is larger than a critical threshold, grains 
are dislodged and carried away so that the porosity of the medium is altered by erosion. This 
in turn affects the local hydraulic conductivity and pressure in the medium and results in the 
growth and development of channels that preferentially conduct the flow. Our multiphase model 
involves a dynamical porosity field that evolves along with the volume fraction of the mobile and 
immobile grains in response to fluid flow that couples the spatiotemporal dynamics of the three 
phases. Numerical solutions of the resulting initial boundary value problem show how channels 
form in porous media and highlights how heterogeneity in the erosion threshold dictates the form 
of the patterns and thus the ability to control them. 



The dynamics of fluid flow through porous continua is 
relevant over many orders of magnitude in length scale 
with applications that range from large scale flow through 
fractured rock in aquifers and oil reservoirs to small scale 
flows in natural and artificially engineered tissues and 
gels [IJE]- I n a ^ these cases, flows are characterized by 
large variations in hydraulic conductivity of the medium. 
This heterogeneity is usually ascribed to processes as- 
sociated with the process of consolidation of the porous 
medium via the agglomeration of grains (in geology) and 
cells (in biology), but may also arise due to channels 
that develop in frangible porous structures as fluid flows 
through them. Indeed flow-induced erosive processes on 
the surface of porous media have been implicated in the 
formation of patterns on planetary [3 , littoral [4] and lab- 
oratory [5] scales that involve both unconsolidated and 
consolidated media [6]. 

However these erosive instabilities can also occur in 
the bulk of fluid saturated materials where they can lead 
to internal channelization via the dynamic coupling of 
flow and hydraulic conductivity. Here we address the 
dynamical evolution of channels via flow-induced erosion 
within a saturated porous medium. Our theory for the 
active co-evolution of the phases in a porous medium 
is qualitatively different from the single phase diffusive 
models for the evolution of free surfaces by aggregation 
and erosion [7] or multiphase bulk theories for multiple 
fluids and/or elastic solids interacting with each other 
[8] but yet combines features of both in considering the 
fluid-induced erosion and deposition processes acting in 
the bulk of a solid skeleton. 

To enable a continuum field description of the process, 
we consider a representative volume much larger than 
the grain/pore size with <j) s (x, y, z, t) the volume fraction 
of the immobile solid phase, (j) g {x,y, z,t), the volume 
fraction of the granular mobile phase, and (f)i(x,y : z,t) 
, the liquid volume fraction in the medium so that 
<t>s + <\>g + <\>i = 1- When the flow- induced stress in a 
frangible porous medium exceeds a local failure thresh- 



old, particles are dislodged and mobilized pQ. This leads 
to a local increase in the hydraulic conductivity and a de- 
crease in the local fluid stress, even as the eroded material 
is carried away. Simultaneously, deposition of mobile par- 
ticles can act to decrease the porosity and reroute fluid 
flow. Thus, for a given flow rate, the hydraulic conduc- 
tivity evolves in space and time as a function of both the 
flow-induced stresses and the initial heterogeneity in the 
porosity of the medium, and can lead to complex ero- 
sional and deposit ional patterns. Volume conservation 
for the individual phases (each of which are assumed to 
be incompressible) implies that 

dt^s = -e + d (1) 
dt<t> 9 = +e-d-V-(^u g ) (2) 
dtfo = -dt{<l>s + <t> 9 ) = -V • (^m) (3) 

where e is the rate of erosion of the immobile phase, d the 
rate of deposition of the mobile phase, and u g , ui are the 
velocities of the granular and liquid phases, respectively. 
We note that adding equations (l)-(3) yields the global 
continuity equation 

V • (0^u g + f ui) = 0. (4) 

We will assume that u g = ui = u, i.e. the granular 
and liquid phases have the same velocity and the effects 
of inertia and body forces associated with sedimentation 
are negligible, a reasonable approximation for slow flows 
of nearly jammed grains. Then the continuity equation 
reduces to V • </>u = 0, where <\> = 0^ + 0/. Erosion of 
<p s can occur only when the fluid-induced stress exceeds 
a critical threshold a = cr(0 s ), where <p s = V~ x j v (j) s dv. 
This form of the threshold function characterizes the non- 
local nature of elastic stress distribution and failure in 
the porous medium in terms of a regional average of (j) s . 
Then we may write the local erosion rate e as 

e = k e ^ s ((7- 1 Vp) 2 - a((/> a )) > 0, (5) 

where k e is a rate, and 7 a nominal pressure gradient 
(based on the fluid flow rate and the hydraulic conduc- 
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FIG. 1: Numerical solution to equations (l)-(7) with a steady 
specific discharge q = qo prescribed at the lower boundary 
y = 0, and constant pressure at the upper boundary y = 0.32. 
(A) Spatial distribution of porosity <fi at t = 1 and (B) at t=6. 
(C) (7 _1 VP) 2 (red) and a (black) are plotted along a-a (at 
y = 0.16) at t = 1 and (D) at t = 6 corresponding to the 
panels above. Erosion occurs where (7 _1 VP) 2 > a. At early 
times, this occurs at several locations. As erosion progresses, 
the pressure gradient drops, heterogeneity in a increases, and 
erosion is limited only to the channels. (E) The flux in the y- 
direction, v(<fi g + (pi) plotted at section a-a at t = 0.5 (black), 
t = 1.5 (red), and at t = 6 (blue). (F) As the flux increases in 
eroded regions, it decreases in non-channelized regions. Here 
we use the tank profile for the erosion threshold (see text), 
and n, = n 2 = 1. 



tivity). The form of the erosion rate follows from consid- 
erations of symmetry: a hydrostatic pressure p cannot 
lead to erosion, but a gradient in p can. However, the 
sign of the gradient is not important, so that we have 
chosen the simplest analytic dependence consistent with 
this symmetry [11]. We assume a = 0.5(tanh(27r((/> s — 
0.6)) + 1),0 < a < 1 to mimic the sharp dependence of 
the failure stress on the volume fraction, although later 
we will consider other forms as well. A simple model for 
the local deposition rate the rate at which the mobile 
granular material is converted back to the immobile solid 
phase, is given by 
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FIG. 2: Result of varying the model parameters on the final 
distribution of porosity <j> plotted at t = 6. The scaled specific 
discharge q specified at y = is (A) decreased to 0.3 and 
(B) increased to 2 Large values of q lead to complete erosion, 
whereas weak q leads to no erosion. (C) The erosion rate k e is 
increased by a factor of 10, i.e. Hi =0.1, and (D) deposition 
rate kd is increased by a factor of 100, i.e. Hi = 100. 



The form of the deposition with a rate kd is based on 
a binary collision picture - mobile grains will come to 
rest only if they interact with immobile grains, with a 
threshold <jf s . In the porous medium, we assume that 
the volumetric flow rate per unit cross-sectional area, i.e. 
the specific discharge q = u(0/ is given by Darcy's 

law[T2] 
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q = ucj) = —D((j))Vp, where D = 



(t> 6 l 



(7) 



d = k, 



I) 
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where d > 0. 



(6) 



Here, D ((/)), the hydraulic conductivity, is assumed to fol- 
low the Carman-Kozeny relation pQ and is in general a 
nonlinear function of the local fluid (pore) volume frac- 
tion (f) = <pi + (j) g . Here l g is the nominal pore size (which 
scales with the grain diameter), \i is the dynamic vis- 
cosity of the interstitial fluid. Using l g ~ 1mm and the 
dimensionless constant A ~ 10 2 for spherical grains [1 
results in D ~ 10~ 5 m 3 s/kg in a water saturated medium, 
which we assume to be isotropic. Erosion leads to a wide 
range in 0, and consequently, variations in D. 

We use the domain size, L = lm, specific discharge 
qo = lcm/ 's, and time T = L/qo [9 to make the prob- 
lem dimensionless, so that the dimensionless parameters 
in the problem include the thresholds for erosion and 
deposition a and 0*, as well as the ratio of deposition 
to erosion rates III = and the ratio of advection to 
erosion rate n2 = -^3. We solve equations (1-7) numeri- 
cally in 2-dimensions (x, y) using a finite volume method, 
with a constant scaled specific discharge q = (0, q) at the 
inflow boundary y = 0, while at the outflow boundary 
y = L y we set pressure p = (atmospheric pressure). 
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FIG. 3: The final distribution of porosity (at t = 6) is shown 
to be sensitive to the initial heterogeneity in a arising from 
(j) s . The upper row shows the porosity distribution resulting 
when the initial distribution of <j> s is set to 0.8 everywhere, 
except in specific places where it is decreased by 1% from the 
uniform background value (A) along two lines, and (B) along 
10 lines, each being one grid cell wide. In the lower row, the 
standard deviation (sd) in the initial perturbation to cj) s is 
varied from its previous value of sd — 0.01. (C) sd — 0.001. 
(D) sd = 0.03. 



In the x direction, we use periodic boundary conditions 
at x — and x = L x , with a square domain of dimen- 
sion L x = 0.32, L y = 0.32 and a uniform grid resolution 
A = 0.05. Prescribing the inlet pressure instead of the in- 
let discharge leads to either no erosion (if (7 -1 Vp) 2 < a ) 
or catastrophic erosion if the pressure gradients are larger 
than the threshold. Thus we prescribe a specific dis- 
charge at the inlet, and a vanishing pressure at the out- 
let. The pressure is determined by iteratively solving the 
Poisson equation V(-D(0)Vp) = obtained by substitut- 
ing (7) into (4), then calculating the erosion rate e and 
the deposition rate d and finally evolving equations (1)- 
(3) to update the volume fraction of the phases <fr 3 , (j) g , <pi 
from one time step to the next, with time step At = 10 -4 . 
Starting with an initial mean volume fraction of mobile 
grains <j) g = throughout the domain and a mean liquid 
volume fraction <\)\ = 0.2 with an additive white Gaussian 

noise (standard deviation sd = \/ (</> 2 ) — (<fri) 2 =0.01) 
which leads to variations in a in the medium, we allow the 
system to evolve until it reaches a quasi-steady state. The 
non-local form of the erosion threshold cr(0 s ), where the 
overbar denotes a weighted spatial average is calculated 
numerically as (j) s (i,j) = 0.2(j) s (i,j) + 0.1(j) s (i ± 1, j =b 1) 
(thus averaging over the 8 surrounding neighbors of a 
grid point) and amounts to a radius of influence of the 
fluid-induced stress that extends a few grain diameters. 



In Fig. 1A,B we show two snapshots of a porous 
medium being evolved with q = at the inlet bound- 
ary, and III = II2 = 1, as it starts to erode inhomoge- 
neously and reaches its final channelized state. Dynam- 
ically, this process arises via positive feedback: flow is 
enhanced through the regions of low solid fraction (high 
hydraulic conductivity) as the fluid scours out a channel, 
while regions with a higher solid fraction (and strength) 
and lower hydraulic conductivity are circumvented by 
the flow. Indeed, in Fig. 1C,D we see the interplay be- 
tween the heterogeneity in the erosion threshold a and 
the squared pressure gradient (7 -1 Vp) 2 , leading to an en- 
hanced flow through regions of high hydraulic conductiv- 
ity at the expense of low flow through other regions with 
the passage of time (Fig. IE) even as the total flow rate 
remains the same. In Fig. IF, we provide a global view 
of the process: erosion continues until the average poros- 
ity over the entire domain increases sufficiently so that 
the pressure gradient everywhere falls below the thresh- 
old for erosion and the system approaches a quasi-steady 
state, where erosion and deposition become vanishingly 
small everywhere. 

To understand how this steady state depends on the 
dimensionless parameters, we first vary the scaled specific 
discharge q at the inlet. When q < q c , i.e. II2 < II C , a 
critical scaled flow rate that depends on the initial poros- 
ity distribution in the medium, no erosion or channeliza- 
tion occurs, because the pressure gradients everywhere 
are smaller than the erosion threshold. For II2 = II C , a 
single narrow channel and secondary partial channel are 
formed as shown in Fig. 2A. As II2 is increased further, 
the number of channels as well as the width of channels 
increases (Fig. IB); for even higher flow rates, the entire 
medium begins to erode away as shown in Fig. 2B. In all 
cases, the mean steady-state porosity increases with the 
specific discharge, linearly at first, before it asymptotes 
slowly to a steady state. Varying the erosion and deposi- 
tion rates also leads to variations in the erosion patterns; 
in Fig. 2C, we show the erosional pattern resulting from 
a 10- fold increase in k e (III = 0.1), which leads to the 
faster evolution of channels. Conversely, increasing kd 
100-fold so that III = 100 increases deposition and lead- 
ing to the blockage of channels (Fig. 2D), although the 
average number or size of channels does not change in 
either case relative to when IIi = 1 (corresponding to 
Fig. IB). 

A question that naturally arises is the mechanism for 
the selection of channel spacing and width when fluid 
flows through a nominally homogeneous porous medium. 
The natural length scales in the problem are the system 
size L the nominal pore size l p which evolves with time, 
but remains a microscopic length, and the length scales 
Qo/k e ,qo/kd; the latter control the dynamical evolution 
of the channels but not their final state. Increasing the 
domain size does not affect the number or size of chan- 
nels. What remains is the threshold for erosion a; since 
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FIG. 4: The evolution of the porosity is sensitive to the func- 
tional form of the erosion threshold a. Here the final distribu- 
tion of porosity cj) is shown for other choices of a. (A) a = <p a , 
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(B) a — <j> 3 . These results can be compared with Fig. 1, where 
a = O.5(tanh(27r(0 s — 0.6)) + 1) (where subtracting 0.6 instead 
of 0.5 provides a slight asymmetry to the tanh profile.) 

the onset of channelization is strongly influenced by fluc- 
tuations in the porosity (and thus the fragility) of the 
medium, we expect that linear stability analysis of the 
base state should predict that channels form at locations 
where a is smallest initially. Thus for the same inlet 
specific discharge, the size and number of channels is a 
function of a(x, y, 0). In Fig. 3A we show that for a given 
inlet specific discharge, if <j(x, 0) = f(x) has a single 
minimum, a single channel forms and grows until the 
pressure gradient falls below the erosion threshold, while 
in Fig. 3B, we see that if a(x : y,0) has multiple minima, 
multiple channels form. Of course, it is not sufficient 
to consider the mean value of the threshold; instead one 
must account for the complete probability distribution of 
the erosion threshold. For our simple Gaussian model of 
disorder, if the variance in the threshold for erosion (or 
equivalently the porosity fluctuations) is also changed, 
this leads to variations in the patterns as well. In Fig. 
3C,D, we show that an increase in the standard deviation 
of the initial white noise leads to greater heterogeneity 
in the channel number and spacing. 

Finally, we consider the functional form of the erosion 
threshold cr((j) s ). In Fig. 4 A we see that for a = <p s the 
final morphology of the erosion patterns is more uniform 
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compared to that shown in Fig. 4B for a = </> s , which 
is itself less variable than for the case when a follows 
a tanh profile (Fig. 1-3). We thus see that the form of 
the erosion threshold function, and its initial, possibly 
heterogeneous, spatial distribution are crucial in deter- 
mining the growth and form of the channels, which a 
simple linear analysis alone cannot capture. 

Our theory of flow-induced channelization in porous 
media focuses on the simplest facets of the phenom- 
ena that involve changes in porosity, pressure gradients 
and flow in a minimal multiphase theory that involves 
fluid, granular and immobile phases interacting with each 
other. Below a critical flow rate, little or no erosion oc- 
curs. Above this threshold, the porous medium starts 
to erode heterogeneously at locations where the critical 



threshold is lowest; positive feedback then enhances ero- 
sion locally in other frangible regions leading to oriented 
regions of higher porosity that are the hallmarks of chan- 
nels. 

Recent experiments with bidisperse granular mixtures 
in a Hele-Shaw cell [9] are qualitatively consistent with 
the strong dependence of the erosional patterns on ini- 
tial heterogeneities in the strength of the porous medium. 
A natural next step is to compare experiments with 
theory quantitatively, and in particular to understand 
the physics associated with the microscopic parameters 
a, k e , kd in laboratory and terrestrial physical systems 
and possible generalizations to biological systems in the 
context of vascularization in natural and artificial set- 
tings [TO] . 
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